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A biological regulatory network can be modeled as a discrete function / that contains all available 
information on network component interactions. From / we can derive a graph representation of the 
network structure as well as of the dynamics of the system. In this paper we introduce a method to 
identify modules of the network that allow us to construct the behavior of / from the dynamics of 
the modules. Here, it proves useful to distinguish between dynamical and structural modules, and 
to define network modules combining aspects of both. As a key concept we establish the notion 
of symbolic steady state, which basically represents a set of states where the behavior of / is in 
some sense predictable, and which gives rise to suitable network modules. We apply the method to a 
regulatory network involved in T helper cell differentiation. 

1 Introduction 

Qualitative methods present a rigorous mathematical framework for modeling biological systems for 
which experimental data needed to determine kinetic parameters and mechanisms is lacking. The com- 
ponents of the system are modeled as variables adopting only finitely many values, so-called activity 
levels. In the simplest case, we obtain a Boolean representation, where the values and 1 may for ex- 
ample represent a gene being inactive or active. In the general case, each component can have several 
activity levels, which may be appropriate depending on the biological data, and often is useful when 
modeling components that influence several other network components. A vector assigning each com- 
ponent an activity level then represents a state of the system. The information about network structure as 
well as the logical rules governing the behavior of the system in state space is represented by a discrete 
function /. 

Although discrete networks are a strongly simplified representation of the original system, complex 
networks are hard to analyze, not least because the state space grows exponentially with the number of 
components. So, methods to reduce the complexity of the analysis are of great interest. One approach 
is to deconstruct the network in smaller building blocks that can be analyzed more easily, which leads to 
the notion of network modularity. 

The idea of decomposing networks into modules is well-established in systems biology, although 
the notion of network module is not clear-cut. Often modules are defined based on biological criteria, 
that have to be translated into mathematical properties in order to identify them in a mathematical model 
(see e.g. [5, 4, 12]). Other approaches focus purely on the the graph representation of the network 
structure. Modules are defined as subgraphs satisfying graph theoretical characteristics often related to 
connectivity [3], or with statistical significance in comparison with random networks [1 1, I]. In addition 
to this structural view, there are also approaches to find dynamical modules, see e. g. [6], that focus 
on identifying behavioral characteristics. However, in general the results obtained by analyzing such 
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modules in isolation do not translate to the original network, since additional influences have to be taken 
into account once the module is re-embedded in the original system. Here, the key is finding conditions 
that allow to draw conclusions about a complex network from knowledge obtained from module analysis, 
as e. g. possible in the modular response analysis approach in the context of metabolic networks and 
steady state fluxes [9, 2]. 

In this paper, we focus on the discrete modeling approach, presenting a method to identify network 
modules that allow us to derive precise information on the dynamics of the original system from the 
results of the analysis of the modules, building on ideas and significantly extending results from [17, 16]. 
In particular, we show that we can explicitly construct attractors of the original systems from network 
module attractors. Here, modularity is a key concept, and we exploit a purely structural as well as a 
purely dynamical view of modularity to eventually determine network modules combining important 
aspects of both. The core notion in our method is that of symbolic steady state. Such a state represents 
a set of constraints on the activity levels of the network components that allows us on the one hand to 
focus on dynamics restricted to subsets of state space, on the other hand enables us to identify dynamical 
and structural modules that render the basis for defining suitable network modules. 

The paper is organized as follows. In the next section we describe the discrete modeling formalism 
we use throughout the paper, and introduce structural, dynamical and network modules. In Sect. 3 we 
establish the notion of symbolic steady state as well as related concepts. This is followed by the main 
results concerning network analysis utilizing modules in Sect. 4. We then illustrate the results for a class 
of networks, namely networks with input layer. In Sect. 6 we apply the method to the analysis of a 
regulator network involved in T helper cell differentiation proposed in [10]. We close with conclusions 
and perspectives. 

2 Discrete regulatory networks 

In this paper we model regulatory systems as discrete functions which capture all available information 
about network interactions and the logical rules governing the behavior of the system. Throughout the 
text, let us consider a system consisting of n € N network components oii, . . . , a„. In the following 
we identify a component a, with its index i to simplify notation. Each component is interpreted as a 
variable which takes integer values that represent the different activity levels of the component. The 
literal meaning of those levels may be very different for different network components, for example they 
can represent levels of substance concentration, gene activity, presence or absence of a signal and so on. 
A vector assigning each component an activity level represents a state of the system, and the dynamics 
of the system is represented by state changes due to component interactions. 

Definition 2.1 For all i G {1, . . . ,n}, let pi € N, and set = {0, 1 , . . . , pi}. Set X = Xi x ■ ■ ■ x X„, and 
let f = {fi, . . . ,fn) : X ^ X be a function. We call f a network comprising n components. For each 
j G {1, . . . the value pt is the maximal activity level o/a,, andXi is called the range of at. The setX 
is called the state space off. 

Each coordinate function /,■ of / describes the rules governing the behavior of the i-th network component 
depending on the state of the system. But / carries not only dynamical but also structural information 
on the system. Both aspects can be represented by directed graphs derived from / as we will see in the 
following two subsections. In the remainder of the paper / denotes a network as introduced in Def. 2. 1 
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2.1 Structure 

We represent the structure of a network by a signed directed (multi-)graph, where vertices represent the 
network components, and an edge from a,- to aj signifies that the value of fj depends on the activity level 
of a,-. The sign of the edge represents the character, i. e., activating or inhibiting, of the interaction. This 
description is inherently local in nature, so we first introduce a structural representation depending on the 
state of the system. This notion was introduced for Boolean functions in [13] and is used for multi- value 
functions in the form considered here in [14]. 

Definition 2.2 Let x ^ X. By G{f){x) we denote the directed signed (multi- )graph with vertex set V = 
{«!,...,«„} and edge set E{x) C.V xV x {+,—}. An edge {i,j,£) belongs to E{x) iff there exists 
c,- G {— 1 , + 1 } such that x,- + c; G Xj and 



where the function sgn : Z {+, — ,0} satisfies sgn(O) = 0, sgn(z) = + if z > 0, and sgn(z) = — if 
z<0. We call G{f) {x) the local interaction graph of / in x. 

The local interaction graph in x is closely related to the discrete Jacobian matrix as introduced in [15] in 
the Boolean case. Note that in the multi-value other than in the Boolean case it is possible that G{f) (x) 
contains parallel edges. There are at most two parallel edges from one vertex to another which then have 
opposite sign. 

The local definition is easily extended, if we are interested in a representation of the interactions 
influencing the system behavior in larger subsets of state space. 

Definition 2.3 Let Y C.X. We denote by G{f){Y) the union of the graphs G{f){x), x We denote the 
graph G{f) {X) also by G{f) and call it the global interaction graph of /. 

In Fig. 1 (b) we see the global interaction graph of the network defined in Fig. 1 (a). The local interaction 
graph in the state (1,1,0) is shown in Fig. 2 (a). To simplify notation we often write G{x) and G{Y) 
instead of G(/)(x) and G(/)(7), respectively, if the corresponding function / is clear from the context. 

When analyzing interaction graphs we are in particular interested in modules of the graph, a term 
for which there exists a variety of definitions as mentioned in the introduction. For our purposes it is 
convenient to use the term in the broadest sense, initially. 

Definition 2.4 A directed (multi-)graph G' = {V^ ) is called a subgraph of a directed (multi-)graph 
G = {V^,E^) ifV^' C ^G' ^ g^g^ gjgg £G' f^^g fj^ff^ end-vertices in V^'. We call a 

subgraph ofG{f) structural module of /. 

Note that for a structural module the vertex set, the edge set or both may be smaller than for G(/). 
For example, the graph shown in Fig. 1 (d) is a structural module of the function / given in the same 
figure. In general, local interaction graphs are structural modules of /. Of course, this definition is not 
very useful for analyzing the network structure or finding characteristics of the system. However, it is 
of conceptual advantage for us in the endeavor of defining network modules that combine structural and 
dynamical characteristics. 
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(a) 

/ : {0,l,2}x{0,l}x{0,l} -> {0,l,2}x{0,l}x{0,l} 

/l(.x) = 0. if ,Xl = .T2 = A .T3 = 1, 

fi{x) = 2, if XI = 2 A X2 = 1, 

fi{x) = 1, otherwise, 

h{x) = 1, if (Xi > 1) V (X2 = 1), 

f2{x) = 0, otherwise, 

/3(a;) = X2 A ^X3 



Figure 1 : A network (a), its global interaction graph (b), and its state transition graph (c), where attractors 
are indicated by fat gray edges. Logical disjunction, conjunction and negation are represented by V, A, 
and -I , respectively. In (d) a structural, in (d) a dynamical module of /. 

2.2 Dynamics 

There are different approaches to deriving the dynamics of /. Commonly used is the so-called syn- 
chronous update strategy, where the successor of a state x is its image under /. A more realistic assump- 
tion is that not all changes indicated by differences in component values of x and f{x) take the same 
amount of time to be executed, since they may represent very different biological processes. However, 
we lack the information to decide which of those processes of activity level change is the fastest. There- 
fore, all possible state transitions are taken into account leading to a non-deterministic representation of 
the dynamical behavior. Furthermore, we assume that a component value changes only by absolute value 
one in each transition, even though the function value may indicate a bigger change. This update method 
is called asynchronous update [18, 19]. 

Definition 2.5 We denote by S{f) the directed graph with vertex set X and edge set E(S{f)) defined 
as follows. An edge {x,x') is in E{S{f)) for states x = {xi,. . . ,x„),x' = {x[,. . . ,xJJ G X if and only if 
xf = f{x) =xorx'j= Xi + sgn(/,(x) —Xi)for some i G {1, . . . ,n} satisfying Xi ^ fi{x), andx'j = xjfor all 
j ^ i. We call S{f) the asynchronous state transition graph of /. 

To analyze state transition graphs we use, in addition to standard terminology from graph theory such as 
paths and cycles, the following concepts. 

Definition 2.6 An infinite path (x^\x\ . . .) in S{f) is called trajectory. A nonempty set of states D is 
called trap set if every trajectory starting in D never leaves D. A trap set A is called attractor if for all 
x^ ,x^ G A there is a path from x^ to in S{f). Attractors of cardinality greater than one are called cyclic 
attractors. A state x is called steady state, if there exists an edge x — ^ x, i. e. if f{x) = x. 

It is easy to see that each trap set contains at least one attractor, and that attractors are the terminal 
strongly connected components of S{f). They represent asymptotically stable behavior and often have 
clear biological meaning. 
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In Fig. 1 (c) we see the state transition graph for the network introduced in (a). The system has two 
cyclic attractors, namely {(1, 1,0), (1, 1, 1)} and {(2, 1,0), (2,1, 1)}. 

As for the structural graph, we can define modules of the state transition graph as subgraphs, i. e. sets 
of states and corresponding state transitions representing fractions of the system's dynamics. However, 
it may also be of interest to only focus on the behavior of a subset of network components, which we can 
derive from the state transition graph by projection. 

Definition 2.7 Let S' = {Y,E{S')) be a subgraph ofS{f). Let -.X ^ Iliei^i projection on the 

components in the ordered .ve? 7 C {1, . . . ,«}. We define Tl\S') as the graph with vertex set Tl^{Y) and 
edges n\v^) k\v^) for v\v^ G 7 such that there exists an edge — > in S{f), and = or 
= + sgn(/,(v') — vl)for some ieL We call n'{S') dynamical module of /. 

Any subgraph of S{f ) is a dynamical module of /, whith I in the above definition chosen as the set 
{!,..., n}. Fig. 1 (e) shows a dynamical module of the function / given in the same figure. Here, we 
choose the subgraphs' consisting of the cyclic attractor {(1,1,0), (1,1,1)} and the corresponding edges 
between the two attractor states. Then the dynamical module 7i^^'^\s') of / is the graph in Fig. 1 (e). 

Again, we do not incorporate any restrictions in the definition that ensure a significance of the mod- 
ules, as e. g. in the approach in [6] where the authors focus on projected dynamics that are independent 
of the behavior of the rest of the system. 

2.3 Network modules 

As mentioned in the introduction, we are looking for subnetworks of / that are on the one hand easier 
to analyze than / itself, and on the other hand carry information of importance for understanding the 
original system. We define modules of the network / as follows. 

Definition 2.8 We call a fiinction g a network module of /, if there exist Y C. X and an ordered set 
7 C {1, . . . ,n} such that g : n^Y) — >■ 7t'{X) satisfies 

^ °f\Y = goTi'lr, 
where |y denotes the restriction of a function to the set Y. 

We call g autonomous, if there exist integer intervals Zi = {a,- , a, + 1 , • • • , fe;}, ai < bi, for a// i G { 1 , . . . , k}, 
k = card7, such that n\Y) = Zi x • • • x and ifg{n'{Y)) C n'{Y). 

Let us again illustrate the notion using the example introduced in Fig. 1. For 7 = {(1,1,0), (1,1,1)}, 
we have /(y) = {(1, 1,0), (1, 1, 1)}. If we set 7 = {1,3} and ^ : {1} x {0,1} ^ {0,1,2} x {0,1} with 
g((l, 1)) = (1,0) and g((l,0)) = (1, 1), then n'{f{x)) = g{n'{x)) for x G Y. Since g{n'{Y)) C n'{Y) 
holds, g is an autonomous network module of /. Note that the set 7 is a set on which the behavior 
of / can in some sense be characterized by the behavior of the components in 7. If we add e. g. the 
state (1,0,0) to the set Y in our example, then there is no function g satisfying o f\Y = go;r'|y, 
since ;r'(l,0,0) = 7r^(l,l,0), but 7r^(/(l,0,0)) = (1,0) / (1,1) = 7r^(/(l, 1,0)). That is, we cannot 
distinguish the behavior of / in states (1,0,0) and ( 1 , 1 , 0) if we only have information on the components 
in 7. 

In general, network modules represent rather local aspects of the network in the sense that they 
describe the influences acting on a subset of components in a set of states. However, the information 
inherent in a network module does not necessarily suffice for determining dynamics beyond a single 
transition step. The second condition for autonomous network modules g allows to apply g iteratively 
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on states in 7t'{Y), while the first ensures that we can derive trajectories according to the asynchronous 
update rule in a projection of state space. Moreover, the first condition allows to apply Def. 2.2 to g. Thus, 
for autonomous network modules g we can determine an interaction graph G{g) and a state transition 
graph S{g). By abuse of notation we also denote G{g) the graph derived from G{g) by renaming the 
vertices I, . . . ,k o{ G{g) with the indices in / while preserving the order. This allows us to identify the 
interaction graph of g with a subgraph of G(/) . 

Lemma 2.9 Let g be an autonomous network module as introduced in Def. 2.8. Then G{g) is a structural 
and S{g) is a dynamical module off. 

Proof As already mentioned, we associate each / G / with /' G { 1 , . . . , ^} via an order-preserving map- 
ping, and rename each vertex of G{g) with indices in I according to this mapping. Obviously, the vertex 
set of G{g) resp. S{g) is a subset resp. a projection via 7i' of a subset of the vertex sets of G(/) resp. 
S{f). Let / G /, and choose /' G {1, . . . as above, i. e., n' maps the i-th component of a state x G X to 
the Z'-th component in 7l\X). Then we have /,(>') = {7i\f{y)))ii = gii{n\y)) for all y G F. Application 
of this equation to the conditions defining edges in Def. 2.2 and 2.5 easily renders that each edge in G{g) 
is also an edge of G{f), and that S{g) = 7i\S'), where 5' denotes the subgraph of S{f) with vertex set Y 
and edges y^ y^ of 5(/) with G F. □ 

For the network module g as defined as illustration for Def. 2.8 the state transition graph is shown in 
Fig. 1 (e). We rename the vertex set {1,2} of G(^) with / = {1,3}. The graph G(g) then consists of the 
vertices a\ and as and a negative loop on as shown in Fig. 1(d). 

In the following sections, we focus on developing a method to determine network modules useful in 
the analysis of /. 

3 Symbolic steady states and frozen components 

Often network components are involved in a number of specific tasks. Thus, although a network com- 
ponent may have a large range, only subsets of the range may be of interest when focusing on specific 
network behavior. To exploit this observation, we introduce the following notation. Here, we call the 
set [a,-, ft,] := {a,-, a,- + 1, . . . — 1,Zj,} C X; a discrete interval, if a,- < hu with [a,, a/] := {a;}. In the 
following, we identify {a,} with a, for all «; G X;, / G {1, . . . ,«}, and call a, regular value. We will use 
intervals of cardinality greater than one instead of regular component values, if we do not have enough 
information to determine the exact component value. Following the terminology in [17, 16], we call 
intervals [a;,^/] with a, < bi symbolic values. 

We now need to integrate symbolic values in the dynamical analysis. Here, we generaUze ideas from 
[17, 16]. 

Definition 3.1 For every i <E {1, . . . ,n} letXP denote the set {[ai,bi] C X,- | a,- < bi} of discrete intervals 

in the range X,. Set = x • • • x X^. We call elements in X regular, elements in X^ \ X symbolic 
states. By J{M) we denote the set of all symbolic valued components ofMforM G X^. Define 

F -.X^ ^X^,M^iFiiM),...,FJM)) with FAM) = [mmfi{x),m&xfiix)] for all jG{l,...,n}. 
We call a state M G X'-' \X satisfying F{M) = M symbolic steady state. 
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The elements of X are subsets of X. The functions / and F coincide on the set X of regular states 
which we identify with the elements of of cardinality one. In general, if a component function value 
fi{M) is regular, this means there is enough information inherent in M to exactly specify its value, while 
a symbolic value represents the fact that we have not enough information to do so. However, it may be 
possible to at least derive some constraint for the function value represented by the interval boundaries. 
For our running example given in Fig. 1 the state ([1,2], 1, [0, 1]) = ([1,2], l,Xj,) = F(([l,2], IjX^)) is a 
symbolic steady state, where we can determine F2(([l,2], 1, [0, 1])) = 1 exactly, obtain the constraint that 
the first component cannot have value 0, but have no information on the third component. 

We are particularly interested in regular components of a symbolic state M that remain fixed on all 
trajectories starting in M. Keeping in mind that we consider the asynchronous update strategy, we can 
find a superset of the set of states reachable from M by the following procedure. We define M° := M 
andM^ := [min(Mj-i UF,(M*-i)), max(Mj-i UFj{M''-^))] for all keN. Since the boundaries of the 
intervals Mj decrease resp. increase monotonously and are bounded by resp. the maximal activity level 
Pj, the sequence {M^)keNo converges to a symbolic state M representing a superset of the set of from 
M in S{f) reachable states. In particular, no trajectory starting in M can leave M. We call M extended 
forward orbit ofM. The next definition is in reference to the notion of frozen cores in random Boolean 
networks introduced by S. Kaufman [8]. 

Definition 3.2 Let / G { 1 , . . . , «}. IfM is a symbolic state with regular component Mi such that Mi = M,- 
for the extended forward orbit M ofM, then we say that the i-th network component is a frozen component 
of M, or a component frozen to value M,-. The set I of all frozen components of a symbolic state M is 
called frozen core ofM, and is denoted by {I,M). 

If a component j of a symbolic state M has symbolic value Xj, then of course the j-th component of the 
extended forward orbit is also Xj. For the example network in Fig. 1 the symbolic state M = {Xi, 1,X3) 
coincides with its extended forward orbit. Thus, the frozen core of M is given by ({2}, (Zi, 1,Z3)). 

Clearly, the frozen core of a symbolic steady state coincides with its set of regular components. 
Moreover, we can use the frozen core of a symbolic state to obtain a symbolic steady state, as the next 
statement shows. 

Theorem 3.3 Let {I,M') be the frozen core of a symbolic state M' G X^. Set M^ := M' and M^ := 

F{M'^^^) for all ^ G N. Then (M^)yteN converges to a regular or a symbolic steady state M. We call M 
the (symbolic) steady state derived from {I,M'). 

Proof If the sequence converges to a limit M, clearly F{M) = M follows from the definition of the 

sequence. The state M is a regular or a symbolic steady state depending on the cardinality of M. We 
show convergence of {M'')^^^ by proving via induction that the sequence is decreasing monotonously 
with respect to the subset relation. That is, we show M'"*"^ C m' for all / G N. 

For i G /, we have Mj = Mf by the definition of frozen components. If i G J{M^) = {1, . . . ,n} \ /, we 
have Mj C Mf by the definition of the extended forward orbit. 

Now, let Z G N and assume C M''-^ for all k<l. Recall that M/+i = F(M'), = [min^g^//,(x), 
max^gji^//,(x)] andM? =F{M'^^)i = [min^g^(-i /,(.x:),max^g^/_i /,(;c)] forall/G {1,. ..,«}. Since C 
M'"^ we have min^g^/-i fi{x) < min^g^/ fi{x) < max^^^i fi{x) < max^g^/-i fi{x), and thus M.+^ C M? 
forall jG {l,...,n}. □ 

As mentioned above, for our running example ({2}, {Xi,l,X^)) is the frozen core of the state {Xi,l,Xi), 
where Xi = {0, 1,2} mdX^ = {0, 1}. Since the state coincides with its extended forward orbit, we start 
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Figure 2: Consider / as given in Fig. 1. In (a) the local interaction graph G((1,1,0)) of /. In (b) 
the graph G(/|m) for the symbolic steady state M = ([1,2], 1, [1,0]). Here, heavier gray edges indicate 
the two components of G^(M). In (c), (d) and (e) the state transition graphs S{f^'), S{f^^) and S^, 
respectively. 

the sequence with = (Zi, 1,^3) = ([0,2], 1, [0, 1]), and = ([1,2], 1, [0, 1]) = M is the symbolic 
steady state derived from ({2}, {Xi,l,Xj)). 

4 Network analysis using modules 

In the following we want to determine network modules, such that the results of the module analysis can 
be directly used to obtain information on the behavior of the original system. Clearly, identification of 
such modules is generally only possible if we exploit at least some coarse knowledge of structural as 
well as dynamical characteristics of the original system. It turns out that the information inherent in a 
symboUc steady state M is sufficient to determine network modules. We proceed by first associating a 
dynamical module with M, then derive a structural module, and finally define a network module suitable 
for utilization in network analysis. 

Let us start by analyzing dynamical characteristics associated with a symbolic steady state M. If 
X G M is a regular state, then f{x) G M by definition of F. More precisely, xt + sgn{fi{x) — xi) G M,-, since 
the interval bounded by Xi and fi{x) is a subset of M,-. Thus, every trajectory starting in M remains in M. 
We formulate this fact in the following statement. 

Proposition 4.1 IfM is a symbolic steady state, then the set of regular states represented by M is a trap 
set. 

By definition, the subgraph of 5(/) with vertex set M and all edges in 5(/) with both end-vertices in 
M is a dynamical module. The result above shows that this module is of significance for the dynamical 
analysis of /, since a trap set always contains at least one attractor. 

In order to associate a structural module with M, we have to recognize some important properties 
attached to the regular components of M. While symbohc components may be dynamically active in the 
trap set M, i. e., the components change their values along at least some trajectories in M, the regular 
components remain fixed regardless of the behavior of the symbolic components in M. This means that 
the network components with symbolic values, i. e., components in J{M), do not influence the behavior 
of the network components with regular values in the trap set M. In turn, the influence of the regular on 
the symbolic components remains the same for all states in M. This motivates the following definition 
describing the structural modules associated with M. 



Definition 4.2 Let M be a symbolic steady state. By G^{M) we denote the (multi-)graph with vertex set 
V^{M) = J{M) and edge set E^{M) = {{i,j,e) G E'^ \ i,j G J{M)}, where denotes the edge set of 
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G(/|m)- We call a graph Z = {V^,E^) component of G^{M), if the undirected graph derived from Z is 
a maximal connected subgraph of the undirected graph derived from (M). 

Note that we use the global interaction graph instead of the local interaction graph G(M) in 

the definition. The difference is that we only consider edges derived from component value changes 
in M instead of in X (compare Def. 2.2 and 2.3), and thus capture all interactions functional in M. In 
particular, there are no edges originating in frozen components of M, since it is not possible to vary their 
value without leaving M. 

To illustrate the above notions let us again consider the example introduced in Fig. 1. We have seen 
that the state M = ([1,2], 1, [1,0]) is a symbolic steady state. In Fig. 1 (c) we can easily see that the set 
of regular states x G M is a trap set which contains both attractors of the system. The global interaction 
graph G(/jM) is shown in Fig. 2 (b). We obtain G^(M) simply by eliminating the vertex Ui- The two 
components of G^{M) are the loops originating in a\ resp. 0:3. 

The components of G^ (M) are the structural modules we associate with M. In preparation for our 
definition of a network module derived from a structural module we need to verify that network com- 
ponents belonging to different components of G^ (M) do not influence each others behavior in M. This 
property is captured in the following lemma, which has already been proved under shghtly different 
conditions in [17, 16]. 

Lemma 4.3 Let M be a symbolic steady state, and letZ\,...^Zk be the components ofG^ (M). Consider 
a union Z of arbitrary components Zj. Let x,y such that x,- = yifor all i ^ Z. Then fi{x) = fiy) 
for all i ^ Z. In particular, for M' G such that M- = Mi for all i ^ Z and M[ C Mi for i G Z, we have 
Fi{M') = Fi{M) = Mi = M'i for all i ^ Z. 

Proof For / ^ J{M), i. e., for a frozen component / of M, we have fi{x) = fi{y) = Fi{M) since x^y G M. 

Let i G J{M) \ Z, and assume fi{x) / fi{y)- We know that if xj 7^ yj then j G Z. Since x,y EM, we 
can define a sequence {x = x^ ,x^ , . . . = y) in M such that x' and x'+^ differ in one component only, 
which is in Z, and the corresponding component values differ by absolute value one. Since /,(x) / fi{y), 
it follows that /,(xO / /,(;c'+') for some / < m. According to Def. 2.2 there exists an edge in G(/|m) 
from some component in Z to i, which is a contradiction. □ 

The lemma shows that the frozen core of M constitutes a boundary between the components of G^ (M) 
that enables us to analyze their behavior in isolation from each other. To do so, we now derive a network 
module from a component Z of G^ (M) by defining a function 

Lemma 4.4 Let M be a symbolic steady state, and let Z = {V^, E^) be a component ofG^ (M). Set = 
cardy^. Let be an order preserving bijectionfrom {l,...,k}to V^. SetX^ = M(Z(i) x • • • x M^z^j^y We 
define f^:X^^ X^, f^ = 7i^oFo p^, where p^:X^^ X° with pf{z) = Mi for i ^ Z and pf{z) = ZiZ(,) 
for i G Z, and : X^ X^ is the projection on the components ofZ. 

Then f^ is an autonomous network module, and is called the network module of / derived from Z. 

Proof The function maps regular states to regular states, since Z is disjoint from J{M) \Z in G^{M) 
and thus Fi{p^{z)) G X,- for z G X^, i G Z, according to Lemma 4.3 and the definition of F. We now show 
n^o f\M = f^on^\M. Let j ^ {!,... ,k} andx G M. To simplify notation we drop the superscript Z from 
;r^, and l^. We have FiQ)(p(;r(x))) = ^i(;)(x) = fi{j){x) according to Lemma 4.3 and the definition 
ofp andF. Thus, (;r^o/|M(x))y = /,(,.)(x) =F,(,)(x) =F,(,.)(p(;r(x))) = (;roFop(;r(x)))y = /f(;r(x)). 
Furthermore, since M is a symbolic steady state X^ = n{M) satisfies the conditions regarding the domain 
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of an autonomous network module given in Def. 2.8. □ 

It is easy to see that the global interaction graph G(/^) is isomorphic to Z, so the structural module 
derived from M matches the one derived from The dynamical modules derived from all functions 
Z component of G^{M), i. e., the state transition graphs S{f^), constitute a breakdown of the coarse 
dynamical module, i. e., the trap set M, which we used to determine first the structural and then the 
network modules. We end this section by showing that these finer dynamical modules are building 
blocks of the dynamics of /. More specific, we show that we can compose the state transition graph of / 
from the state transition graphs of the network modules and the frozen components of M. Again, in 
the following we generalize results from [17]. First, we define the composition of the graphs S{f^). 

Definition 4.5 Let M G he a symbolic steady state, and let Z\,..., Zf„ be the components of (M). 
We then denote by the graph with vertex set M and edge set E^. An edge — > x^ belongs to the edge 
set iff 

x^ =x^, and n^' (x^) ^ Tt^' belongs to S{f^i ) for all j G { 1 , . . . , m} , 
or if there exists j E {1 , • • . ,m} such that 

n^^{x^) Tt^'{x^) is an edge in S{f^-') and x] = xj for all i ^ V^J . 
We call the product state transition graph corresponding to M. 

The next theorem confirms that the method of composing the state transition graphs of the network 
modules renders the subgraph of S{f) derived from the state set M. 

Theorem 4.6 Let M G be a symbolic steady state, and let Zi, . . . be the components of G^(M). 
Let S{f) \m denote the subgraph ofS{f) with vertex set M and all edges in S{f) with both end-vertices in 
M. Then S{f)\M = S'^ . 

Proof Recall that M is a trap set, and that for G M there is an edge x — > ^ in S{f) if and only if 
x' = f{x) =x orx- = x; + sgn(^(x) — x;) for some i G {1,. . . ,n} satisfying x,- ^ fi{x), andx^- = Xj for all 

First, we note that fi{x) = x, = M, for all x G M and i ^ 7(M). Furthermore, for i G /(M) there 
exists k' G {l,...,m} such that i G V^*^', and there exists /' such that = i, with l^v being the 

bijection introduced in Def. 4.4. As seen in the proof of Lemma 2.9, we have fi{x) = /^*' (tt^*' (x)). 
Therefore, x G M is a fixed point of / iff 7r^^(x) is a fixed point of f'^i for all j G {1, . . . ,m}, and 
Xi + sgn(/,-(x) —Xi) = {n^ki {x))ii + sgn(/^*' {n^ki {x)) — {n^k' (x))/,). Thus, we can construct each edge in 
S{f ) \m from an edge in some S{f^') and vice versa. It follows from Def. 4.5 that |m = 5^- D 

The reasoning in the proof of Theorem 4.6 leads immediately to the following statement. 

Corollary 4.7 Let M and Zi, . . . ,Zm be as in Theorem 4.6. For all i & {\,. . . , m} let A,- be an attractor 
in S{f^'). Then A := {a G M | Vj G {1, . . • ,m} : Tl^^{a) G Aj] is an attractor in S{f). Moreover, every 
attractor in S{f)\M can be represented in this manner as product ofattractors in S{fzj), j G {1, . . . ,m}, 
and component values Mi for i ^ J{M). 

We illustrate the results on our running example from Fig . 1 . The graph (M) f or M = ( [ 1 , 2] , 1 , [0 , 1 ] ) 
has two components, Zi consisting of a positive loop on ai, and Z2 being a negative loop on 0:3, as can 
be seen in Fig. 2 (b). We define /^i : [1,2] ^ [1,2] as in Def. 4.4 with tt^^ ((MJ ,M^,M^)) = M[ for all 
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(a) 

/:{0,1,}3^{0,1}3, 

fi{x) = Xl, 
f2{x) = Xi, 
fsix) = {xi A X2) V -iSs 



(b) 



a2 



aiJ 



+ 



(c) (0,1,0) 
(0,1,1) 



or) 



(1,1,0) 

A 



(1,1,1) 



(0,0,0) 



(0,0,1) 



(1,0,0) 
(1,0,1) 



Figure 3: A network with input layer in (a), its global interaction graph in (b), and its state transition 
graph in (c). 



M' CM and p^' (z) = (z, 1, [0, 1]) for all z G [1,2]. This reduces the function fi given in Fig. 1 (a) to 
f^^{z) = z. The state transition graph S{f^^) is given in Fig. 2 (c) and consists of two steady states. 
Analogously, we obtain S{f^^), which comprises a single attractor of cardinality 2 as shown in Fig. 2 
(d). Applying Def. 4.5 we obtain the graph which is shown in Fig. 2 (e). Comparison with the state 
transition graph S{f ) given in Fig. 1 (c) illustrates Theorem 4.6 and its corollary. 



5 Networks with input layer 

In this section we introduce a class of networks, for which we can easily find a set of symbolic steady 
states such that all attractors of the original system can be constructed from the attractors of the network 
modules derived from the symbolic steady states. We extend results obtained in [17]. 

Definition 5.1 We call f a network with input layer, if there exists i G {1, . . . ,n} such that ft = idx;. A 
vertex satisfying this condition is called input vertex. 

For every input vertex a,-, we have no incoming edges except a positive loop in the global interaction 
graph. However, this structural criterion is not sufficient for identifying an input vertex. Networks with 
input layer are well-suited for modeling, for example, signal transduction networks, where receptors may 
be modeled as input vertices. 

In the following we assume that / is a network with input layer. Without loss of generality we assume 
that ,o;<:, ^ e { 1 ,...,«}, are the input vertices of /. 

We can immediately note one important property of /. The frozen core of a symbolic state M' G 
with M- e Xi for i E {I,. . . ,k} and M- = Xi for i > k is given by ({1, • . .,k},M'), since M' is its own 
extended orbit. In other words, we can easily derive a (symbolic) steady state from each combination 
of values for the input vertices using Theorem 3.3. The resulting set of symbolic and possibly regular 
steady states is sufficient to determine all attractors of /. 

Theorem 5.2 Let A be an attractor of f. Then there exist input values x, G X,-, i E {1, . . . ,A:} such that 
either A=M is a regular steady state, or we can construct A from M as shown in Cor 4. 7, where M is the 
fixed point derived from ({1,. . . ,^},M'), M' G X^ withM[ =Xifor j G {1, . . . andM[ =Xifor i > k. 

Proof Since ai . . . a^: are input vertices, we have a, = a ■ for all i G { 1 , . . . , A:}, and we set M- = a, for the 
input vertices. Next, we show that for every state xeX with x, = M- for all ? G { 1 , . . . , k}, there exists a 
trajectory leading to M. 
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Figure 4: Global interaction graph of the Th cell differentiation network introduced in [10]. Arrows 
represent activation, crossed lines represent inhibition. 

Let (M');eNo with := M' and M' := F{M'^^) for all Z G N be the sequence converging to M 
introduced in Theorem 3.3. Recall that M' C M'"' for all / G N. Let £M'\M. Then there exists 
I G No such thatx^ G M' \M'+'. Then, according to the definition of F and since (M')/gpjQ is decreasing 
with respect to the subset relation, either all successors of in S{f) axe in M' \M'+^ or there exists a 
successor in M'"'"^ In the latter case, we label that successor a:^ . Otherwise, we repeat the procedure for 
all successor of x^, check again and if necessary repeat again. Since all images of all states in M' lie in 
M'+' and since the state space is finite, we eventually find a state G M'+' such that there exists a path 
fromx*^ tox^ with all states of the path except x^ lying inM' \M^~^^. Since (M')/eNo is converging to M 
we can thus construct a path from each state in M' to M. 

It follows that there is no trap set, and thus no attractor, in M' \M. If M is a regular steady state, then 
M = A is the only attractor in M'. Otherwise A is a composition of attractors of network modules derived 
from M and the frozen components of M as shown in Cor 4.7. □ 

We illustrate the results on the simple Boolean network given in Fig. 3 which has one input vertex, namely 
«! . We start calculating the two (symbolic) steady states from the input values xi = and xi = 1 . In the 
first case, we get F((0, [0, 1], [0, 1])) = (0,0, [0, 1]) and F((0,0, [0, 1])) = (0,0, [0, 1]), i. e., (0,0, [0, 1]) is 
a symbolic steady state, and the associated structural module is a negative loop originating in a3. When 
looking at the corresponding network module, we obtain the function : [0, 1] [0, 1], f^{z) = ^z, and 
the derived attractor for the original network is the set { (0, 0, 0) , (0, 0, 1 ) } . For the input value xi = 1 we 
get the sequence ((1, [0, 1], [0, 1]), (1, 1, [0, 1]), (1, 1, 1), (1, 1, 1), . . .), so the procedure renders a regular 
steady state of the system. As can be seen in Fig. 3, the regular steady state and the attractor derived 
from (0,0, [0, 1]) are the only attractors of the system. 

6 Analyzing Th cell differentiation 

T helper cells, short Th cells, are important players in the vertebrate immune system. They can be 
sub-classified in Thl and Th2 cells, which are involved in different immune responses. Both originate 
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Table 1 : Coordinate functions and ranges for the components of the Th cell network. 



from a common precursor, promote their own differentiation and inhibit proliferation of each other. In 
[10] L. Mendoza proposes a model for a control network of Th cell differentiation consisting of 17 
components, 13 of which are represented by Boolean variables while the remaining 4 components have 
three activity levels. The logical rules governing the behavior of the system are given in Table 1 and the 
global interaction graph can be seen in Fig. 4. Note that the model depicted in Fig. 4 and Table 1 differs 
slightly from the model introduced in [10], namely we altered the logical functions associated with the 
vertices lFN-j3, lL-12 and 11-18. In Mendoza's model, the three components are modeled with constant 
functions with value 0, representing the wild type in some sense (see [10]). The constant values are 
changed when considering specific artificial environmental conditions. We model these vertices as input 
vertices and consider for the wild type the situation where all input values are zero. Clearly, the attractors 
of both models coincide. Modeling lFN-j3, lL-12 and lL-18 as input vertices also makes sense from a 
biological point of view since all three vertices represent substances not reproduced by Th cells. If we 
want to mimic experimental conditions where cells are cultured in media saturated with one or more of 
these substances, we can easily do so by focussing on the part of state space where one or more of the 
input vertices have value one. Note that we have no further input vertices in our model. 
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For the wild type, i. e. the situation where all input values are set to zero, Mendoza identifies four 
attractors all of which are fixed points of the function / given in Table 1 . Each one has a clear biological 
interpretation [10]. We now want to apply our analysis technique using symbolic steady states to the 
wild type. 

We fix the values of the input vertices to zero and as a first step determine the corresponding symbolic 

steady state, that is, the symbolic steady state derived from the frozen component set ({1,2,3}, xi =X2 = 
X3 = 0). Iterating the state M0:= (0,0,0, [0,1], [0,2], [0,1], [0,2], [0,1], [0,1], [0,1], [0,1], [0,2], [0,1], [0,1], 
[0,1], [0,2], [0,1]) we get 

Ml := /(M«) = (0, 0, 0, 0, [0, 2] , [0, 1] , [0, 2] , [0, 1] , 0, 0, [0,1], [0, 2] , [0, 1] , [0, 1] , [0, 1] , [0, 2] , [0, 1] ), 
m2 := /(Ml) = (0,0,0,0, [0,2], [0, 1], [0,2], [0, 1],0,0, [0, 1], [0,2], [0, 1],0,0, [0,2], [0, 1]), 

/(M2)= M2. 

We obtain a symbolic steady state with 8 regular components, and no further constraints on the re- 
maining components. The local interaction graph G^(M) is shown in Fig. 5 (a). Analysis of the cor- 
responding subnetwork renders four fixed points, namely (0,0,0,0,0,0,0,0,0), (1,0,1,0,1,1,0,1,0), 
(2,0,1,0,1,1,0,2,0), (0. 1.0,1,0,0,1,0,1) GX5 xXexX? xXg xXii XZ12XZ13 xXi^xXn. The steady 
states of the original network derived from these fixed points match the four steady states found in [10]. 
The state space of the original model consists of 663552 states. Fixing the input values still leaves us 
with 82944 states to consider. The state space of the structural module associated with the symbolic 
steady state M^ contains only 2592 states. 

Not all of the combinations of input values render a significant simplification of the network analysis. 
In the worst case, for example if we choose input values xi = 0, X2 = X3 = 1, we can only derive the value 
for X4 but no further constraints on structure and behavior of the system. This in itself is of course an 
interesting observation from a biological point of view, since in that case we can deduce that cross- 
regulation plays an important role at an early stage of signal transduction. 

On the other hand, some combinations of input values lead to very small network modules. Let us as a 
last example consider the input values x=l,X2=X3 = 0, representing an overabundance of IFN-j3 . Start- 
ingwithMO:= (1,0,0, [0,1], [0,2], [0,1], [0,2], [0,1], [0,1], [0,1], [0,1], [0,2], [0,1], [0,1], [0,1], [0,2], [0,1]) 



we get 






Ml 




/(M«) = 


(1,0,0, 1, [0,2], [0,1], [0,2], [0,1],0,0, [0,1], [0,2], [0,1], [0,1], [0,1], [0,2], [0,1]), 


M^ 




/(Mi) = 


(1,0,0, 1, [0,2], [0,1], [0,2], [0,1],0,0, [0,1], [1,2], [0,1],0,0, [0,2], [0,1]), 


M^ 




/(m2) = 


(1, 0,0,1, [0,2], 0, [0,2], [0,1], 0,0,1, [1,2], [0,1], 0,0, [0,2], [0,1]), 


M^ 




/(m3) = 


(1,0,0, 1, [0,2], 0, [0, 1], 0,0,0, 1, [1,2], [0, 1],0,0, [0,2], [0, 1]), 


M^ 




/(M^) = 


(1,0,0,1, [0,2], 0,[0,1], 0,0,0,1, 1,0,0,0,[0,2], [0,1]), 


M^ 




/(M5) = 


(1,0,0,1, [0,2], 0,[0,1], 0,0,0,1, 1,0,0,0,[0,2],0), 


M'^ 




/(m6) = 


(1,0,0,1, [0,2], 0,[0,1], 0,0,0,1, 1,0,0,0,[1, 2], 0), 


m8 




/(M^) = 


(1, 0,0,1, [1, 2], 0,[0,1], 0,0,0,1, 1,0,0,0,[1, 2], 0), 


m9 




m') = 


(1,0,0,1, [1,2], 0,1, 0,0,0,1, 1,0,0,0,[1, 2], 0), fiM^) = M^. 



We obtain a two-component module consisting of a positive loop on T-bet and an activating edge from 
T-bet to IFN-7 which can be seen in Fig. 5 (b). Both components of the module originally have three 
activity levels, but are both constrained to levels 1 and 2 in the module dynamics. Thus, we only have to 
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Figure 5: Subnetworks of the Th cell differentiation network associated with the symbolic fixed points 
derived from the input values = X2 = JC3 = in (a) and = 1, ^2 = = in (b). 



analyze a state space of cardinality four instead of a state space consisting of 82944 states. The module 
has two steady states, namely (x5,xi6) = (1,1) and (x5,xi6) = (2,2), which translate to two steady states 
in the original network. Again, this is in agreement with the results in [10] (supplementary material). 

Application of our analysis method to this model thus offers two advantages. First, the complexity 
of the analysis of the dynamics is reduced, since we only have to focus on the smaller network modules. 
Furthermore, identification of the modules themselves is of interest, since they represent the part of the 
system responsible for the decision of the system's fate, i. e., which attractor is reached. An interesting 
next step would then be to check whether the mathematically derived network modules coincide with 
subsystems of known biological importance. 



7 Conclusion 

In this paper we introduced a method to analyze discrete regulatory networks utilizing suitable network 
modules. We used the notion of symbolic steady state, which allows us to determine such network 
modules using coarse dynamical and subsequently structural modules derived from the symbolic steady 
state. Lastly, we can associate network modules with the structural modules exploiting the properties of 
the frozen core of the symbolic steady state. We then can construct the dynamics of the original network, 
and in particular its attractors, in a subset of state space explicitly from the state transition graphs of the 
network modules. This paper not only gives a rigorous definition of different aspects of modularity but 
notably extends results in [17, 16]. In particular, the detailed analysis of the Th cell network becomes 
possible because of the refined notion of symbolic steady state. 

A variety of aspects provide possibilities for fruitful future work. Firstly, we want to focus on further 
options for easily computing symbolic steady states. Here, we introduced a method suited for networks 
with input-layer, but the resulting symbolic steady states might not be minimal with respect to the subset 
relation on symbolic steady states, and thus the resulting network modules can possibly be further refined. 
We anticipate results when we focus on certain classes of functions / describing networks, in particular 
(nested) canalyzing functions [8, 7]. A different direction of interest is to use the network modules not 
only for obtaining the system's dynamics but also for a refined stability analysis. Perturbations resulting 
in changes in the dynamics of / might not be noticeable in every network module (see [6]). Such consid- 
erations are of similar interest when using the synchronous update strategy, which is also often utilized 
in discrete modeling, so a translation of our results to synchronous update networks seems worthwhile. 
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Furthermore, we need to compare our results to other well-established modularization techniques that 
aim at reducing the analysis complexity such as the modular response analysis of biochemical networks 
[9, 2]. Although the underlying modeling approaches are very different, some ideas may be transferable 
or render complementing results. Lastly, we plan to apply our methods to further biological examples. 
Here, a comparison of the network modules and their associated structural and dynamical modules with 
subsystems of known biological importance could lead to interesting insights. 
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